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Abstract 


The  performance  bounds  of  a  passive  acoustic  array  operating  in  a 
turbulent  medium  with  fluctuations  described  by  a  von  Karman  spectrum 
are  investigated.  This  treatment  considers  a  single,  monochromatic, 
spherical-wave  source  and  a  line-of-sight  propagation  path.  The 
Cramer-Rao  lower  bounds  of  the  wave-front  angles  of  arrival  are 
calculated  for  an  unknown  parameter  set  which  includes  the  propagation 
distance,  turbulence  parameters,  source  phase,  and  signal-to-noise  ratio. 
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1.  Introduction 


The  Cramer-Rao  lower  bounds  (CRLBs)  of  the  angle-of-arrival  (AO A) 
estimates  for  a  plane  wave  incident  on  a  passive  sensor  array  were 
investigated  in  Refs.  [1,2],  herein  refered  to  as  the  PW  Papers.  In  the  PW 
Papers  it  was  assumed  that  the  wave  had  propagated  through 
atmospheric  turbulence  with  fluctuations  described  by  a  von  Karman 
spectrum.  The  received  signal  was  modeled  as  a  complex  Gaussian 
random  variable  with  a  non-zero  mean.  For  simplicity,  a  single 
monochromatic  source  and  a  line-of-sight  propagation  path  were 
assumed.  The  technique  was  based  on  that  of  Wilson  [3],  which  followed 
the  general  frame  work  of  Song  and  Ritcey  [4], 

The  present  paper  is  a  logical  continuation  of  Refs.  [1,2],  but  for  an 
incident  spherical  wave.  The  theoretical  model  of  Refs.  [1,2]  is  used,  with 
the  appropriate  changes  made  for  a  spherical  wave.  The  CRLBs  of  the 
azimuthal  and  elevational  AOAs  are  calculated.  The  propagation  distance, 
signal-to-noise  ratio  (SNR),  turbulence  parameters,  and  phase  of  the 
source  are  also  treated  as  unknown  parameters.  As  the  estimates  of  the 
AOAs  will  degrade  when  they  are  simultaneously  estimated  with  other 
parameters,  the  couplings  between  the  estimates  of  the  AOAs  and  the 
estimates  of  the  other  parameters  are  also  calculated. 

The  paper  is  ordered  as  follows:  The  theoretical  model  is  discussed  in 
section  2  and  CRLBs  are  formulated  and  discussed  in  sections  3  and  4.  The 
numerical  results  and  concluding  remarks  are  given  in  sections  5  and  6, 
respectively.  Tables  of  the  acronyms  and  symbols  are  given  in  appendix  A. 
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2.  Theoretical  Model 


The  theoretical  model  for  the  probability  density  function  (PDF)  of  a 
spherical  wave  that  has  propagated  through  atmospheric  turbulence  is 
outlined  in  this  section.  The  following  notation  shall  be  used  throughout 
this  paper:  [•]*  denotes  the  complex  conjugate,  [  ]T  the  transpose,  [•]*  the 
Hermitian  adjoint  (complex  conjugate  transpose),  (•)  the  ensemble 
average  or  expectation  value,  and  [  ]^  and  [•]'*  the  vector  components 
parallel  and  transverse  to  wave  propagation,  respectively. 


2.1  Signal  Model 

Consider  an  array  with  N  sensors.  The  total  received  signal  at  each  sensor 
is  taken  to  be  the  sum  of  the  wave  that  has  propagated  from  the  source  of 
interest  and  from  noise.  Both  contributions  are  time  dependent.  The  total 
received  signal  s  may  be  written  as  a  column  vector  with  N  elements,  one 
corresponding  to  each  sensor, 

s(<f>,d,t)  =  p(<f>,6,t) +  n(t)  ,  (1) 

where  p  is  the  wave  that  has  propagated  from  the  source  with  azimuthal 
and  elevational  AOAs  4>  and  0,  respectively,  and  where  n  is  the  noise.  It  is 
assumed  that  the  noise  at  the  sensors  are  mutually  uncorrelated  and  that 
the  noise  has  a  complex  Gaussian  distribution  with  zero  mean  and 
variance  at  each  sensor.  Exact  solutions  for  the  pressure  field  of  the 
source  and  its  PDF  are  not  known,  but  solutions  to  its  moments  are. 
Therefore,  it  is  approximated  that  p  is  complex  Gaussian  distributed  with 
mean  p  and  covariance  matrix  Cp, 

Hi  =  {pi)  and  [Cp]y  =  (pip*)  -  mp* .  (2) 

The  moments  are  determined  from  the  results  in  the  open  literature  for 
acoustic  wave  propagation  in  a  moving  random  medium  as  discussed  in 
the  following  section.  We  assume  that  p  and  n  are  uncorrelated  so  that  the 
total  signal  s  is  Gaussian  distributed  with  mean  p  and  covariance 

C  =  Cp  +  o*IN.  (3) 

This  signal  model,  in  which  the  real  and  imaginary  parts  are  Gaussian 
random  variables  with  equal  variances,  is  reasonable  for  strong  or  weak 
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scattering  in  the  presence  of  strong  diffraction  (the  Rytov  extension 
region).  It  is  less  well  suited  to  situations  where  both  scattering  and 
diffraction  are  weak  (geometric  acoustics),  in  which  case  the  phase 
variance  dominates  the  signal  behavior  [8], 

2.2  First  and  Second  Moments 

The  pressure  field  associated  with  a  sound  wave  propagating  in  a  moving 
random  medium  is  characterized  by  a  dosed  set  of  fluid  dynamic 
equations.  The  small-angle  parabolic  and  Markov7  approximations  may  be 
used  to  obtain  the  statistical  moments  of  the  sound  field  in  closed  form. 
These  approximations  are  valid  in  far  field,  for  small  scattering  angles, 
and  for  £  »  A  >  i,  where  A  is  the  wavelength,  and  £  and  £  are  the  outer 
(integral)  and  inner  (Kolmogorov)  length  scales  of  the  turbulence, 
respectively.  We  use  the  solutions  for  the  first  and  second  moments  of  the 
pressure  field  as  given  by  Ostashev  [5],  who  generalized  the  results  in 
Refs.  [6,7]  to  include  fluctuations  in  the  medium  velocity.  The  solution  for 
the  second  moment  is,  however,  valid  for  normal  incidence  across  two 
sensors. 


2.2.1  Normal  Incidence 


Consider  a  sound  wave  that  is  propagating  with  wave  number  k  =  ker, 
fc  =  27r/A,  where  A  is  the  wavelength  of  the  source.  Let  the  observation 

point  be  n  =  [z,  yu  Zi]T,  so  that  x  »  =  (y?  +  zf)1/2.  From  Ref.  [5],  the 

first  moment  at  i\  for  an  incident  spherical  wave  is 


Mi  =  Ph  (r i)  e  yx  ,  (4) 

where  7  is  the  extinction  coefficient  for  the  first  moment  and  pH  is  the 
sound  field  in  the  absence  of  random  inhomogeneities 

P»(ri)  =  Pi  ei4>i ,  Pi  =  —  ,  <K:  =  kn  +  7 ,  (5) 

ri 

where  y  is  the  phase  of  the  source,  A  is  the  pressure  at  r,  =  r0  (A  and  r0 
are  real-valued  constants). 

The  small-angle  approximation  (SAA),  [i.e.,  parabolic  approximation],  for 
a  spherical  wave  propagating  in  free  space  is  simply  the  Taylor  series 
expansion  up  to  first  order  in  1  /x 


exp  (iAx;) 

n 


—  exp  (iA:z)  exp 


(6) 
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Consider  now  the  observation  points  r *  =  [x,  y,,  ^]T  and 

Tj  =  [x,  yj,  Zj]T,  so  that  x  »  i?j,  From  Ref.  [5],  the  second  moment  for 

normal  incidence  across  sensors  i  and  j  is 

(p*Pj  >  =  Ph  (r<)  Ph  (rj)  e-o(^)x ,  (7) 

where  p^  =  r,  —  r,-  is  the  sensor  separation  vector  (transverse  to  the 
propagation  direction),  a  is  the  extinction  coefficient  for  the  second 
moment,  and  pH  (ritj)  are  given  by  the  parabolic  approximation 
[equation  (6)].  The  extinction  coefficients  for  the  first  and  second  moments 
are  related  by 

7  =  a  (oo)  /2 .  (8) 


2.2.2  Oblique  Incidence 

We  now  wish  to  derive  approximate  expressions  for  the  first  and  second 
moments  for  oblique  incidence.  Let  the  vector  from  the  center  of  the  array 
to  the  source  be  r  =  rr .  The  azimuthal  and  elevational  AOAs,  <b  and  9,  are 
measured  with  respect  to  the  center  of  the  array,  so  that 
r  =  [cos  (f>  cos  6 ,  sin</>cos0,  sin0]T.  Let  r'  =  [x',  y',  z']T  be  the  vector  from 
the  center  of  the  array  to  the  ith  sensor,  r*  =  r  -  r[  be  the  vector  from  the 
ith  sensor  to  the  source,  and  =  r,  17  =  r'  —  r'  be  the  vector  between 
the  jth  and  ith  sensors.  An  illustration  is  given  in  figure  1.  (When 
comparing  results  between  the  plane-wave  case  and  the  spherical-wave 
case,  keep  in  mind  that  in  this  paper  and  in  Ref.  [2],  6  is  the  elevation, 
whereas  in  Ref.  [1],  9  is  the  declination  and  9'  is  the  elevation.) 

For  oblique  incidence,  we  must  take  care  in  deriving  approximations  for 
the  moments.  A  consistent  treatment  of  both  the  phases  and  attenuations 
of  the  moments  is  necessary  to  ensure  that  the  covariance  matrix  Cp  is 
non-singular.  For  normal  incidence,  the  second  moment  is  calculated 
using  the  SAA,  where  pH  (rz)  is  approximated  by  expanding  about  ry  up  to 
first  order  in  1  /ry .  For  oblique  incidence,  we  expand  pH  (r j)  about  r  up  to 
first  order  in  1  /r.  In  doing  so,  we  assume  that  the  attenuation  of  the  wave 
across  the  array  is  constant  and  may  be  approximated  by  that  at  the  array 
center.  Therefore  for  pH  (r *)  =  Pi  exp  (i$i),  we  approximate 


Pi  = 


Ar0  Ar0  _ 


=  P0 


=  X  +  «  X  +  k 


r  —  r 


r'  +  rA 


>2 


,/\2 


(r  •  r') 


2r 


(9) 

(10) 


For  the  first  moment  at  the  ith  sensor,  we  approximate  that  the  attenuation 
of  the  wave  due  to  the  random  medium  is  also  constant  across  the  array 
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z 


Figure  1.  Coordinate  system.  The  closed  circles  represent  the  sensors  and  the  open  circle 
represents  the  source.  The  azimuth,  <f>t  and  elevation,  9f  are  defined  with  respect  to  r . 


and  may  be  approximated  by  that  at  the  array  center.  Thus 

Mi  «poel#ie_7r ,  (11) 

where  4>t  is  given  by  equation  (10), 

The  second  moment  for  oblique  incidence  we  approximate  to  be 

(PiPj)  ~  ple^^e~a{pii)r  <S>13  =  ,  (12) 

where  again  is  given  by  equation  (10).  We  are  thereby  assuming  that  for 
every  i  and  j, 

a(pij)  ~  a(Pij)  •  (13) 

The  factor  of  r  in  the  attenuation  term  due  to  the  medium  is  necessary  to 
ensure  consistency  between  the  moments.  In  other  words,  as 


[Cplij  =  ( PiPj >  -  (Pi)  < Pj  > 

in  the  limit  pl3  — ►  oo. 


■Po 


[C 


'PHj 


►  Po  [e_27r  -  e_27r]  ei<t,J 


-  e-2^r]  ei$u  , 

(14) 

=  o, 

(15) 

as  it  should.  Because  of  these  approximations,  we  limit  our  investigation 
to  nominally  normal  incidence  at  a  planar  array. 

We  define  the  y-element  of  the  mutual  coherence  function  (MCF)  matrix 
to  be  the  positive  square  root  of 

_  (PiPj)(p?Pj) 


r2 _ 

13  (piPi)(pjP*j) 


(16) 
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For  the  treatment  here. 


r«  =  |<KP‘)!/p§  =  e-“('>»>'.  (17) 

The  minimum  value  of  the  MCF  occuring  for  ptj  =  oo  is  rmjn  =  We 
define  a  vector  of  the  phases 

5=  [ei$1,  e1*2,  •••  ,  ei*"]T  (18) 

and  the  matrix 

S  -  s®s^  (19) 

Sij  =  exp  (i<f>jj)  ,  (20) 

where  ®  is  the  (right)  Kronecker  product.  We  may  then  write 

m  =  por^2n  ^  (2i) 

cp  =  p20res-P20Tmins,  (22) 

where  0  is  the  Hadamard  product  (element-by-element  multiplication). 
This  mathematical  expression  is  useful  for  computational  purposes.  As  we 
are  approximating  that  the  attenuation  of  the  sound  wave  is  constant 
across  the  array,  we  may  think  of  5  and  S  as  a  steering  vector  and  a 
steering  matrix  that  involve  only  the  phases  of  the  wavefront. 


Turbulence  Model 


The  extinction  coefficients  depend  on  the  structure  of  the  random 
medium.  For  an  incident  spherical  wave  the  extinction  coefficient  for  the 
second  moment  [5,9,10]  is 


a(p)  =  2  nk2 


f 


[/(0)  -  f(pu)]du , 


(23) 


where  /  is  the  two-dimensional  (2D),  or  projected,  correlation  function  for 
the  sound  speed  fluctuations.  For  most  random  media,  including 
turbulence,  a  (p)  initially  increases  monotonically  with  increasing  p,  but 
when  p  exceeds  C,  a  (p)  asymptotically  approaches  a  constant  value.  Since 
/(/>)—>  0  in  the  limit  p  — >  oo,  this  constant  value  is  simply  27,  given  by 

27  =  2nk2f  (0)  =  2 q2k2C ,  (24) 


where  c2  is  the  index  of  refraction  variance.  Hence  the  second  moment 
initially  decreases  with  increasing  p  and  eventually  "saturates"  at  a  fixed 
minimum  value.  Note  that  this  limit  [equation  (24)]  is  the  same  for  both  an 


incident  spherical  wave  [equation  (23)]  and  an  incident  plane  wave 
[equation  (15)  of  Ref.  [1]]. 

The  performance  of  atmospheric  acoustic  sensor  arrays  that  have  a  sensor 
spacing  larger  than  the  height  of  the  array  from  ground  is  affected  by  the 
large  eddies  of  the  energy -containing  (or  source)  subranges  of  the 
turbulence  spectrum.  The  isotropic,  homogeneous  von  Karman  turbulence 
model  describes  the  inertial  subrange  of  the  turbulence  spectrum  more 
realistically  than  the  commonly  used  Gaussian  models,  and  it  still  behaves 
reasonably  in  the  energy-containing  subrange.  The  von  Karman  form  for 
the  2D  correlation  function  is  dependent  upon  the  source  of  the  sound 
speed  fluctuations:  a  scalar  field  is  induced  by  temperature  or  humidity 
fluctuations  and  a  vector  field  is  induced  by  wind  velocity  fluctuations. 
The  2D  correlation  functions  for  a  scalar  field  fs  and  a  vector  field  fv,  may 
be  written  in  the  form  (see  equation  (49)  in  Ref.  [10]  and  equation  (7.112) 
in  Ref.  [5]) 


fs 

fv  (p,  C,0 


2  q2l 

F  (1/3) 

2  q2l 

v^r(i/3) 


(25) 

(26) 


where  l  =  T  (1/3)  £/  [yfnT  (5/6)]  is  a  characteristic  length  scale,  F  ( x )  is  the 
gamma  function,  and  Kv  (x)  is  the  modified  Bessel  function  of  order  v. 

The  MCF  for  an  incident  spherical  wave  is  plotted  in  figure  2  as  a  function 
of  the  index-of-refraction  variance,  q2,  and  the  eharacterstie  length  scale 
normalized  by  the  wavelength,  l/X,  for  both  a  scalar  and  a  vector  von 
Karman  spectrum.  In  presenting  the  results,  it  is  natural  to  use  normalized 
length  scales,  (e.g.,  rj A,  dj A,  etc.),  as  then  the  coherence  has  no  explicit 
wavelength  dependence.  In  (a)  and  (b)  the  MCFs  are  calculated  for 
p/A  =  0.5  and  r/A  =  500.  The  coherence  for  both  spectra  decreases 
significantly  in  the  regions  where  the  index-of-refraction  variance  is  large, 
c2  ~  10-4,  and  the  normalized  characteristic  length  scale  is  small, 

10  <  l/X  <  1.  (The  MCF  is  actually  dependent  upon  the  product  q2r/X- 
therefore,  the  x-axis  alternatively  represents  an  increase  in  r /A  for  fixed 
C2.)  In  (c)  and  (d)  the  same  is  calculated  but  for  p/A  =  3/\/2.  The  larger 
sensor  separation  leads  to  a  more  rapid  decrease  in  the  MCFs  as  functions 
of  the  turbulence  parameters.  For  both  sensor  separations  the  MCF  for  the 
vector  spectrum  is  more  sensitive  to  the  changes  in  the  turbulence 
parameters,  and  its  minimum  with  respect  to  the  turbulence  parameters 
(for  a  fixed  finite  sensor  separation  and  normalized  propagation  distance) 
is  smaller  than  that  for  the  scalar  spectrum. 

In  order  to  compare  the  results  to  those  of  the  PW  Papers,  the  same  has 
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(a)  p/X  =  0.5  (b)  p/X,  =  0.5 


Figure  2.  Coherence  for  an  incident  spherical  wave:  (a)  and  (c)  are  for  a  scalar  von  Kdrmdn 
spectrum;  (b)  and  (d)  are  for  a  vector  von  Karman  spectrum.  All  calculations  are  for  r/\  — 
500. 


been  plotted  in  figure  3  but  for  an  incident  plane  wave.  Note  that  the  MCF 
for  the  spherical  wave  is  larger  than  that  of  the  corresponding  plane  wave. 
Therefore,  based  on  the  results  from  the  PW  Papers,  we  expect  that  the 
CRLBs  of  the  AO  As  for  an  incident  spherical  wave  should  be  smaller  than 
those  for  the  corresponding  plane  wave. 

The  function  rmjn  (the  minimum  value  of  the  MCF  as  a  function  of  sensor 
separation  (p  =  oo)  for  fixed  propagation  distance  and  fixed  turbulence 
parameters)  is  the  same  for  both  a  scalar  and  a  vector  spectrum  and  for 
both  a  plane  wave  and  spherical  wave.  It  is  plotted  in  Fig.  4.  Even  though 
its  value  is  only  dependent  upon  the  product  q2rl/ A2,  it  is  plotted  versus 
the  turbulence  parameters  at  r/A  =  500  for  ease  of  comparison  with 
figures  2  and  3. 
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Figure  3.  Coherence  for  an  incident  plane  wave:  (a)  and  (c)  are  for  a  scalar  von  Karman 
spectrum;  (b)  and  (d)  are  for  a  vector  von  Karman  spectrum.  All  calculations  are  for  rj A  = 
500, 


3.  Formulation 


3.1  CRLB 


The  CRLB  is  calculated  from  the  Fisher  information  (FI),  which  is 
dependent  upon  the  PDF.  Having  completed  the  theoretical  model  of  the 
PDF,  we  now  proceed  to  formulate  the  FI. 


Consider  the  vector  of  real  parameters  ©  =  [0i,  @2,  •  •  • ,  ©a/']T  that  are  to 
be  estimated.  The  FI  for  real  parameters  of  a  complex  Gaussian  probability 
likelihood  function  with  covariance  matrix  C  and  mean  is  [11] 


=  M  tr 


dC  \ 
d©„/ 


+  2  Af» 


/  dtf  t  djx  \ 

Vd©A  dQj  ’ 


(27) 


for  M  independent  and  identically  distributed  data  sets.  If  there  are  N 
sensors  in  the  array,  C  is  an  N  x  N  matrix  and  /x  is  a  column  vector  of 
length  N.  Let  us  use  the  convention  that  A,  u  e  [1,  2,  . . . ,  Af\  are  the 
indices  on  the  parameters  and  i,j  £  [1,  2,  . . . ,  N]  are  the  indices  on  the 
sensors. 

Let  us  define  <r„  =  ^[J-1]^.  We  loosely  refer  to  either  a  or  a2  as  the 
CRLB,  as  the  meaning  should  be  evident  from  the  units  involved.  The 
minimum  value  of  a 2  is  <j10  =  1/  Jw>  (be.,  the  CRLB  when  ©„  is  the  only 
unknown).  As  the  number  of  unknowns  increases,  cr^  will  increase. 

For  example,  suppose  that  there  are  two  unknowns.  For  A  and  u  cyclic 
(e.g.,  if  A  =  1  then  v  =  2), 


where 

Ci2  =  7^-,  0  <  C12  <  1  •  (29) 

Only  if  J12  =  0  does  <y\  =  a\0,  and  the  estimates  of  @1  and  ©2  are  said  to 
be  uncoupled.  As  C12  increases,  a\  increases  from  its  minimum  value  of 
cr^o,  and  a  degradation  of  the  estimates  of  ©1  and  ©2  results.  The  quantity 
£12  thus  provides  a  measure  of  the  strength  of  the  coupling  between,  and 
hence  degradation  of,  the  estimates  of  ©1  and  ©2:  if  C  =  0,  the  estimates 
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are  uncoupled  and  the  CRLBs  retain  their  minimum  values;  if  £12  <C  1,  the 
estimates  of  9\  and  92  are  weakly  coupled  and  the  CRLBs  increase  only 
slightly;  and  if  (12  =  1,  the  estimates  are  fully  coupled,  the  CRLBs  are 
infinite,  and  hence  neither  0i  nor  ©2  can  be  estimated.  It  is  therefore 
advantageous  to  determine  the  conditions  under  which  the  estimates  of 
©i  and  ©2  will  decouple. 

If  there  are  more  than  two  coupled  parameter  estimates,  we  define  the 
coupling  between  the  Ath  and  z^th  parameter  estimates  to  be 


T2 

A  = 

C  JxxJu 


In  this  way,  we  have  a  measure  of  the  coupling  strength  between  any  two 
given  parameters. 


3.2  FI  of  Theoretical  Model 

We  write  the  elements  of  the  mean  vector  as 

K=poe-'V*‘=p„r^e1*<,  (31) 

and  of  the  total  covariance  matrix  as 


where  <3?ij  =  and  is  given  by  equation  (10). 

The  signal-to-noise  ratio  is  related  to  the  noise  variance  by  SNR  =  p^/a2. 

It  is  often  useful  to  express  the  SNR  in  decibels  SNRdB/  SNR  -  iQSNRdB  /1°. 
For  a  spherical  wave  p0  is  dependent  upon  r;  therefore,  we  consider  SNRo, 
the  signal-to-noise  ratio  at  a  distance  72,  as  the  unknown.  Then 
Gn  =  (-4r0)2/ (SNRo7£2).  By  renormalizing  the  FI  by  (_4r0)2,  the  explicit 
value  of  Ar$  is  not  needed. 

The  FI  may  now  be  readily  calculated  from  equation  (27)  for  those 
parameters  we  wish  to  consider  as  unknowns:  <t>,  9,  r,  l,  c2,  and  SNR,.  It 
was  shown  in  the  PW  Papers  that  the  source  phase  y  must  be  treated  as  an 
unknown  parameter  when  a  non-zero  mean  is  considered;  therefore,  y 
must  be  treated  as  an  unknown  for  the  spherical-wave  case  as  well.  For 
brevity,  the  derivatives  of  the  covariance  matrix  and  mean  with  respect  to 
the  unknown  parameter  set  are  not  presented  here. 
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Due  to  scale  discrepancies  between  the  contributions  to  the  FI  matrix, 
numerical  difficulties  may  arise  when  trying  to  invert  J.  As  in  the  PW 
Papers,  we  construct  a  diagonal  matrix  D,  whose  elements  are  the  inverse 
of  the  square-root  of  the  diagonal  elements  of  the  FI  matrix, 

D\„  =  JZ!/2s\w  The  CRLBs  of  the  unknown  parameters  may  then  be 
determined  from 

^ (e„  -  e„)2^  =  [d  (DJD)-1  d]  w  .  (33) 

The  matrix  DJD  may  be  inverted  by  standard  numerical  techniques. 


4.  Discussion 


4.1  No  Turbulence 

Let  us  begin  by  examining  the  case  of  no  turbulence. 

4.1.1  Full  SW 


In  the  absence  of  turbulence,  not  considering  the  SAA  for  the  moment, 

H  =  [pie1*1,  P2e1*2,  pjveI#jV]T  and  C  =  allN  .  (34) 

The  elements  of  the  FI  matrix  are 


Jx 


MN^da^da^  /  2  d$i  d$i  dpi  dPi\ 

d@x  dQu  a*  \Pi  aex  dQv  dox  aeu )  ' 


(35) 


(36) 

(37) 

(38) 

(39) 

(40) 

(41) 
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where 


r  ■  r[  —  x't  cos  6  cos  0  +  y'i  sin  o  cos  0  +  z\  sin  6 
3  (r  * 

-  —  =  —x[  sin  <f>  cos  0  +  y*  cos  <t>  cos  0 


=  —  x\  cos  (j) sin  6  —  ?/•  sin  (j>  sin  0  +  z-  cos  6 . 


For  the  estimates  to  decouple,  the  off-diagonal  elements  must  be  zero. 
Substituting  equations  (42)-(44)  into  equations  (39)-(41)  we  find 


2M  ^  ( A2k2r%r2  A2r%r2\  \(  ,2  ,2\  ^  ^  ^ 

J<t,6  =  — j*  2_^  ( - T-  +  J  [(*<  -!/*  J  cos <f> sin <f> cos 6 sm  6 


+  x\y\  (sin2  (f>  —  cos2  <p)  cos  0  sin  9  —  x'^  sin  6  cos2  0  4-  y\z\  cos  <f>  cos2  6  (45) 

2M  ^  _4.2  fc7'nr  (a?;  sin  d>  cos  0  —  y ,•  cos  d>  cos  0)  .... 

=  -ZT  E -  - 3 - 1 -  <46> 

"  i=  1  i 

2M  ^  A2krlr  (  x\  cos  <j>  sin  9  +  y[  sin  </>  sin  6  —  cos  0)  (  . 

=  -zr  Z. - ^  •  (4/' 

»=i  * 

Thus  in  order  for  the  estimates  of  <f>  and  d  to  decouple  from  the  estimate  of 
X,  for  every  <t>  and  6,  we  must  have 

N  .  N  ,  N  , 

E3  =  53  =  Ej=o-  <«) 

i= 1  *  1=1*  i=i  * 

Therefore,  unlike  the  PW  case,  there  is  no  simple  array  geometry  that  will 
result  in  the  decoupling  of  the  estimates  of  the  AO  As  from  the  estimate  of 
the  phase  angle.  The  conditions  for  the  estimates  of  </>  and  6  to  decouple 
are  also,  in  general,  unattainable  in  practice: 


A  ,  ,  /fc2  1 

Xiyi  [  ,.4  +  r6 


( fc2  J[\  —  V-  ,2  (^L-  4. 

K rf  rf)  ~^Vl  W  ri) 

e*m  (5 
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Suppose  now  that  r  is  also  unknown.  Then 


2M 


' n 


"  i—l  '  i  i  ' 


/\2 


2M  ^  /  A2k27'QT  A2Tq 


Jr<f>  — 


£±v±  / 

^2"  \ 

un  ,-=1  \ 


i=l 

N 


ri 


+ 


°'  1  (r  +  f -r') 


2M  A  / A2fc2rgr  A2r2r 
Jr0  =  -^2 


i—l 


ri 


+ 


rf 


(r  +  r  r') 


A  djr-tj) 

d(f> 

5(f -r') 


d& 


2 M  ^  A2kr^ 


JrX~  a2 


T*V 

n  i=l  A 


(51) 

(52) 

(53) 

(54) 


Because  of  the  factor  r  in  the  second  term  in  parenthesis,  we  see  that  the 
estimate  of  r  will  always  be  coupled  to  the  estimates  of  <j>,  d,  and  y, 
regardless  of  array  geometry. 

Also  note  that  none  of  the  elements  of  the  FI  given  in  equations  (36)-(41) 
or  (51)-(54)  is  dependent  upon  the  value  of  y.  This  is  expected  as  the  value 
of  the  source  phase  should  not  effect  the  estimates  of  the  other  parameters. 
Close  inspection  of  the  second  term  of  equation  (27)  reveals  that  its 
dependence  in  the  FI  should  cancel  regardless  of  whether  we  consider 
turbulence  or  not. 


4.1.2  SAA 


Let  us  now  consider  results  for  the  parabolic  approximation  using  the 
approximations  of  px  and  given  in  equations  (9)  and  (10).  The  elements 
of  the  FI  matrix  are  now  given  simply  by 

T  _  MN  dal  da2  ,  2MN  dpp  dp0  2 Mpg  d$i 

A"  oi  d@x  de„  cr2  OQxdeA  ai  ( 

We  thus  find 


2 Mplk2  A 

cr2  ^ 

n  t=  1 

Jm  =  •£  ' 


J\X  — 


a 

2MNp~ 


n  i= l 

9 


1  + 


1  + 


d(f) 

t:ri\  dJhlil 

86 


(56) 

(57) 

(58) 
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Jr>j‘  - 


a2r2 

un' 


+ 


_  2Mppfc2 


N  r 

E 

i=i  L 

N  r 


1  + 


^  =  ^E[(1+^) 


7  _  tmpQn 

J° *  ~  a2  2^ 

n  i= 1 

_  2Mp%k 2 

_2Mp20k*^ 

Jrd  ~  a2  2^ 

n  i=l 

_  2Mp%k  y^ 
*'rX  "  ,2  L 

n  i=l 


1  + 


2  n  r 

-E 

1=1 

,  ri2  —  ‘  ri)2] 

2r2 

2 

rri 

’(1+ 

f  ■  r\  \  d  (r  • 

r  , 

d4>  J 

r  )  86 

a(f-r')l 

r  ; 

d<f) 

fr'A 

d(f  -r')‘ 

r  j 

ae 

f  r/N 

1  xt 

r  j 

\8(f-  r»' 

)  dcf> 

\-r- 

i - 1 

<M 

"s-T 

CN 

1 

cs 

f  •  r'N 

x  Lt 

r  y 

\  d(r-  r[)' 

/  oe 

\-r- 

f-(f-r')2! 

2r2 

/\2 


1  - 


rf -(»•>{) 

2  r2 


(59) 

(60) 
(61) 
(62) 

(63) 

(64) 

(65) 


Notice  that  the  summations  on  the  lowest  order  terms  are  the  same  as 
those  for  the  PW  case  [compare  to  equations  (43)-(48)  of  Ref.  [1]]. 
Therefore,  we  can  use  the  results  from  the  PW  case  to  minimize  the 
couplings. 

Again  note  that  none  of  the  elements  of  the  FI  given  in  equations  (56)— (65) 
is  dependent  upon  the  value  of  x- 


4.2  Turbulence 


Let  us  now  consider  propagation  through  atmospheric  turbulence. 


4.2.1  Two-Element  Array 

Consider  a  circular  wave  that  is  received  at  a  two-element  array.  Suppose 
that  the  sensors  are  separated  by  d.  In  Appendix  D  of  Ref.  [1],  the  FI  is 
derived  for  a  2-element  array.  The  results  derived  there  are  also  valid  for 
the  consideration  here. 

Let  a  =  C\\,  b  =  \Cn\,  and  c1/2  =  |m|  =  where  pi  and  Qj  are  as 
defined  in  equations  (31)  and  (32).  As  d<&i/d\  =  1  Vi,  from  equation  (D-9) 
of  Ref.  [1],  we  know 
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(66) 


J<fxp  — 


2  Mb2  /d& 


?12 


b2  V  d<t> 


+ 


2  Me 


a 


?h\\ 

a*} 


<9#2 

d<p 


2b 


d<t>  d<j> 


J, 


_  2 Me  (d$ i  |  d§2\ 


^  a  +  b 
_  AMc 

Jxx  ~  JT+6 


d</>  d<t>  J 


(67) 

(68) 


If  d&i/do  =  — 9$2 /d(j>,  the  estimates  of  0  and  \  will  decouple  for  any 
value  of  c.  The  CRLB  of  6  is 


1 


J<t>4>  4X/JXX 


a 


b2 


M  [2 b2  +  c  (a  +  6)]  (9$i2 /d<p) 


2  • 


(69) 

(70) 


J, 


<t>x 


2  Mck 
ci  T  b 

2  Mck 
a  T  b 


Now 


(®i  +  x2)  sin  0  -  (y[  +  y2)  cos  <j> 


(71) 


1  xi~  x2^  Vi  V2  j  •  j.  ,  x'\y>i  t x>2y2  j .  2  j  2  i\ 

^ - — - cos  4>  sin  <t>  -| - (sin  <j>  —  cos  0) 


(72) 


In  order  for  J$x  to  be  zero  for  every  value  of  0  and  c,  the  following  must 
hold 

0th  order  in  1/r  :  x[  =  -x2  and  y[  —  -y2  (73) 

lsf  order  in  1/r  :  x'  *  +  x'2  =  y[2  +  y'2  and  x[y[  =  -x2y2  .  (74) 

The  0th  order  conditions  are  those  found  for  the  PW.  However,  only  one  of 
the  two  1st  order  conditions  in  (74)  can  hold  when  (73)  holds.  Therefore, 
the  estimates  of  <f>  and  \  will  always  be  coupled,  but  this  coupling  may  be 
minimized  by  satisfying  equation  (73), 

If  x\  =  x'2  —  0  and  y[  =  -y'2  —  d/2,  then 
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=  k2d2  cos2  <f> .  (75) 


d$i2 

d<(> 


=  k2 


d (r  •  £1) 
d(j) 


d<t> 


Thus 


T  _  Mc2kd?  cos  4>  sin  <f> 

J<t>x  =  (o  +  6)r 

-  Mr^infcd2sin2</> 

2  (1  +  T  -  2T^n  +  a2/p2)  r 

_ _ M  e~2',r  kd2  sin  26 _ 

2  [l  +  e-a(d/2)r  —  2e-2">T  +  <r2/po]  r 

At  normal  incidence,  4>  =  0,  the  coupling  is  zero.  This  term  decays  by  at 
least  e~2ir /r,  therefore  at  large  values  of  r,  J0X  decays  very  rapidly.  It 
follows  that 

(79) 

(80) 

This  is  the  same  as  found  for  the  PW,  except  that  po  now  has  a  1/r 
dependence.  Note  that  the  lsf  order  corrections  in  the  phase  term 
canceled. 


2  a2  —  b2 

=  W^d?co^^T7{a  +  bj\ 

_  i  ~  r2  ~  2Tmjn  (i  —  rmin)  +  2a 2/pq  (i  ~  rmjn) 

_  Mk2d2  cos2  0  [2 T2  -  3ITmin  +  Tmin  (l  +  a2/p2)  ' ' 


(76) 

(77) 

(78) 


4.2.2  Planar  Array 

For  full  three-dimensional  (3D)  propagation  in  turbulence,  we  assume  that 
the  wave  is  propagating  near  the  x-axis,  take  the  array  plane  to  be  the 
2/2-plane,  and  take  the  origin  to  be  at  the  center  of  the  array.  Numerically, 
we  find:  the  estimates  of  d>,  9,  r,  and  \  are  all  coupled;  the  estimates  <f>  and 
9  are  uncoupled  from  the  estimates  of  t,  c2,  and  SNRo;  and  the  estimates  of 
r,  t,  c2,  and  SNRo  are  all  coupled.  The  results  are  independent  of  the  value 
of  x- 


18 


5.  Results 


The  array  geometry  considered  for  this  analysis  is  a  4  x  4  square  grid  with 
spacing  of  d.  In  all  figures  d/\  =  0.5.  As  the  CRLB  (cr)  for  M  independent 
and  identically  distributed  datasets  is  l/\/M  times  the  CRLB  for  one 
dataset,  all  results  are  presented  for  M  =  1. 

5.1  CRLB  of  AO  As 

As  the  estimates  of  the  AO  As  are  only  coupled  to  the  esimates  of  x,  r/A,  <t>, 
and  9,  only  these  four  parameters  are  considered  as  unknowns  in  this 
section. 

Due  to  the  symmetry  of  the  array,  the  values  of  a 0  and  ae  are  the  same  at 
normal  incidence.  At  4>  =  0,  the  couplings  Qe,  Qx,  and  Qr/X  are  all  zero. 
And  at  6  =  0,  the  couplings  Q0,  Cex,  and  Qr/X  are  all  zero.  Therefore,  at 
4>  =  9  =  0,  =  a$  =  a0o  =  crg0.  All  the  angular  couplings  increase  with 

increasing  <f>  and  9. 

In  figure  5,  a0  for  normal  incidence  is  plotted  versus  ij\  and  c  for 
r/A  =  500  and  SNRo  =  10  dB  at  72 /A  —  500.  A  scalar  von  Karman 
spectrum  is  used.  The  overall  values  of  are  smaller  than  for  the 
plane-wave  case  (see  figure  16  of  Ref.  [1]  or  figure  11  of  Ref.  [2]).  This  is 
expected  as  the  values  of  the  MCF  for  a  plane  wave  are  smaller  than  for  a 
spherical  wave  (refer  to  figures  2-3  on  pg,  9).  The  same  is  plotted  in 
figure  6,  but  for  a  vector  von  Karman  spectrum.  The  ratio  a0/(Ta0,  where 
the  superscripts  refer  to  the  type  of  spectrum,  is  plotted  in  figure  7.  As 
expected,  oQ  for  the  vector  spectrum  is  larger  than  that  for  the  scalar 
spectrum.  As  with  the  PW  case,  the  use  of  a  non-zero  mean  reduces  the 
CRLBs  of  the  AO  As.  For  other  values  of  4>  and  6,  the  behavior  of  CRLBs  of 
the  AO  As  is  similiar,  with  ac,  and  cfq  increasing  with  increasing  o  and  9. 
And  for  other  values  of  the  SNR,  the  behavior  is  similar. 

In  figure  8,  is  plotted  versus  the  normalized  propagation  distance  for 
normal  incidence  and  a  scalar  von  Karman  spectrum.  Two  values  of  c2, 

£/ A,  and  SNRo  evaluated  at  72/A  =  500  are  considered.  In  order  to  see  the 
limiting  behavior  of  the  model,  the  graph  is  extended  to  include  smaller 
values  of  rj A  than  are  valid  for  the  turbulence  model.  At  small  values  of 
r/X,  we  see  that  i0  is  dependent  upon  the  values  of  the  turbulence 
parameters  (particularly  <r)  and  is  independent  of  the  value  of  SNRo  (as  it 
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Figure  5.  CRLB  of  azimuth  as  a  function  of  turbulence  parameters  for  a  scalar  von  Kdrmdn 
spectrum:  normal  incidence,  rf  A  =  500,  and  SNRo  =  10  dB  at  Rq/\  —  500.  Due  to  symme¬ 
try,  cTp  —  do  —  ct^q  —  (Too  at  <j>  —  0  —  0. 


Figure  6.  CRLB  of  azimuth  as  a  function  of  turbulence  parameters  for  a  vector  von  Kdrmdn 
spectrum.  All  other  parameters  same  as  figure  5. 
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Figure  7.  Ratio  of  cr“  to  oA,  as  a  function  of  turbulence  parameters.  Parameters  same  as 
figures  5  and  6. 

should  be  as  SNRo  is  evaluated  at  TZ/X  —  500).  Note  the  difference 
between  the  outward  spherical  wave  propagation  and  the  plane  wave 
propagation  depicted  in  figure  18  of  Ref.  [1]  and  figure  13  of  Ref.  [2]. 

The  same  is  investigated  in  figure  9  but  for  </>  =  0  =  15°.  While  a0  has 
increased  slightly,  the  behavior  is  the  same.  The  coupling  between  4>  and  6 
is  in  figure  10.  For  r/X  =  10,  C©,  e  is  dependent  upon  the  values  of  the  SNR 
and  the  turbulence  parameters;  however,  this  dependence  rapidly 
vanishes  and  C©,  8  becomes  constant  for  r/X  >  100.  Therefore,  for  values  of 
r/X  that  are  consistent  with  the  SAA,  is  independent  of  c2,  l/X,  and 
SNRo.  Figure  11  depicts  the  coupling  between  b  and  \.  For  all  values  of  c2, 
if  X,  and  SNRo,  Cp,  *  diminishes  exponentially  with  increasing  rf  X.  The 
coupling  between  c>  and  r/X  is  plotted  in  figure  12.  Again  for  large  values 
of  ?’/ A,  Cj>,  r/>,  diminishes  exponentially.  The  analogous  graphs  of  figures 
9-12  for  0  (instead  of  o)  have  the  same  behavior. 

The  behavior  of  and  erg,  as  well  as  all  the  couplings,  as  a  function  of 
propagation  distance  is  the  same  for  a  vector  von  Karman  spectrum.  The 
values  of  a0  and  a0  are  slightly  higher.  Refering  back  to  figure  10,  the 
coupling  for  a  vector  spectrum  has  similar  values  for  r/X  <  100  and 
approaches  the  same  constant  for  r/X  >  100. 

The  angular  dependence  of  the  CRLBS  of  the  AO  As  are  shown  in  figure  13 
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Figure  8.  CRLB  of  azimuth  as  a  function  of  normalized  propagation  distance  for  normal 
incidence.  SNR  evaluated  at  71/ X  =  500,  scalar  von  Kdrmdn  spectrum. 


Figure  9.  CRLB  of  azimuth  as  a  function  of  normalized  propagation  distance  for  <j>  =  0  = 
15°.  All  other  parameters  same  as  in  figure  8.  (Analogous  plots  for  ao  have  same  behavior, 
though  smaller  values.) 


Figure  10.  Coupling  between  <j>  and  0  as  a  function  of  normalized  propagation  distance. 
Parameters  same  as  figure  9. 


and  the  coupling  between  the  AO  As  is  shown  in  figure  14  for  a  scalar  von 
Karman  spectrum.  In  13(a-c),  c2  -  1CT4,  i/X  =  10,  r/X  =  500,  and 
SNRo  =  10  dB  at  1Z/ A  =  r/X.  We  see  the  same  symmetry  in  as  in  the 
PW  case,  and  ae  is  again  independent  of  azimuth.  (Refer  to  figures  (15-16) 
in  Ref.  [2]  and  figures  (21-22)  in  Ref.  [1].)  Not  shown  are  Q>,  r/x  and 
which  are  at  most  ~  10-12  at  <f>  =  0  =  15°.  To  within  the  numerical 
accuracy  of  the  calculation,  Q,tX  and  Ce,x  are  both  zero.  In  13(d-f),  all  the 
parameters  are  the  same  except  c2  =  it) -ti.  The  CRLBs  have  decreased,  as 
expected.  Here  C^r/A,  Qj,r/x>  Cp,X’  Ce,x  ~  10~9  at<j>  =  $  =  15°.  The 
coupling  Q  j,  in  figure  14  is  the  same  for  both  cases.  In  fact,  the  coupling 
between  6  and  9  is  the  same  (only  .02  percent  difference)  as  for  the  PW 
case  for  all  parameters  and  both  spectra  when  r/A  ^  100.  The  analogous 
results  for  a  vector  spectrum  have  the  same  behavior. 

It  is  found  for  all  cases,  that  the  use  of  a  non-zero  mean  reduces  the  CRLBs 
of  the  AO  As.  However,  in  the  regions  where  the  AOAs  can  be  estimated, 
this  percent  difference  is  small,  usually  less  than  two  percent.  This  percent 
difference  is  discussed  in  detail  in  Ref.  [12]  for  both  an  incident  spherical 
and  plane  wave. 
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Figure  11.  Coupling  between  <j>  and  \  as  a  function  of  normalized  propagation  distance. 
Parameters  same  as  figure  9.  (Analogous  plots  for  x  have  same  behavior  and  almost 
same  values.) 
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Figure  12.  Coupling  between  <j>  and  r/X  as  a  function  of  normalized  propagation  distance. 
Parameters  same  as  figure  9.  (Analogous  plots  for  Q,  r/\  have  same  behavior  and  almost 
same  values.) 


Figure  14.  Coupling  between  6  and  0.  Results  valid  for  both  spectra  and  all  values  of  c2, 
e/X,  and  SNRo  when  r/X  >  100. 


5.2  CRLBs  of  Other  Parameters 

For  the  PW  case,  it  was  found  [2]  that  there  was  not  sufficient  information 
to  provide  reasonable  estimates  of  the  propagation  distance  and 
turbulence  parameters.  For  the  spherical-wave  model,  we  again  find  that 
we  cannot  estimate  rf  A,  c2,  if A,  or  SNRo-  Even  when  c2,  i/K  and  SNRo  are 
not  considered  as  unknowns,  rf  A  cannot  be  estimated. 

As  the  turbulence  parameters  cannot  be  estimated  from  this  model,  it 
follows  that  they  must  be  calculated  from  models  that  consider 
meteorological  data.  There  are  models  based  on  von  Karman's  spectrum 
that  calculate  the  turbulence  parameters  in  the  energy-containing 
subrange  of  the  turbulence  [9,10].  These  models  consider  contributions  to 
the  sound  speed  variations  from  wind  and  temperature  fluctuations 
produced  by  both  shear  and  buoyancy  instabilities. 
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6.  Conclusions 


We  have  investigated  the  Cramer-Rao  lower  bounds  of  the  wave-front 
angles  of  arrival  for  a  spherical  wave  propagating  through  atmospheric 
turbulence  with  fluctuation  described  by  a  von  Karman  spectrum.  This 
investigation  is  the  logical  continuation  of  previous  investigations  [1,2]  for 
a  plane  wave  source.  Both  investigations  consider  a  deterministic  mean, 
two  bearing  angles,  and  multiple  unknown  parameters.  The  use  of  a 
spherical-wave  model  reduces  the  CRLBs  of  the  AO  As  from  those  of  the 
plane-wave  model  due  to  the  change  in  the  mutual  coherence  function. 
The  coupling  of  the  estimates  of  the  angles  of  arrival  to  the  estimates  of 
the  normalized  propagation  distance  and  source  phase  are  found  to  be 
minimal,  if  not  completely  negligible.  And  the  estimates  of  the  angles  of 
arrival  are  decoupled  from  the  estimates  of  the  turbulence  parameters  and 
signal-to-noise  ratio.  The  coupling  between  the  two  bearing  angles  is  the 
same  as  for  the  plane-wave  model  when  the  normalized  propagation 
distance  is  consistent  with  the  parabolic  approximation. 

As  the  estimation  of  the  turbulence  parameters  does  not  effect  the 
estimates  of  the  AO  As,  it  is  logical  to  extend  the  calculation  of  the  mutual 
coherence  function  to  include  contributions  from  both  scalar  and  vector 
spectra.  Such  a  model  is  considered  in  Ref.  [13],  in  which  the  MCF  is 
calculated  from  meteorological  data. 

The  significance  of  scattering  by  atmospheric  turbulence  is  evident  from 
the  results  presented  here.  In  order  to  understand  and  circumvent 
limitations  on  U.  S.  Army  acoustical  tracking  systems,  it  is  necessary  to 
predict  the  performance  of  acoustic  sensor  arrays  for  various  atmospheric 
conditions.  This  analysis  clearly  demonstrates  the  atmospheric  conditions 
that  are  unfavorable  for  accurate  acoustical  tracking. 

This  investigation  is  limited  by  the  model  of  the  second  moment  for 
oblique  incidence,  as  well  as  the  fact  that  other  physical  phenomena,  such 
as  ground  reflections  and  refraction  by  atmospheric  wind  and 
temperature  gradients,  which  have  not  been  included  in  this  analysis,  may 
have  a  considerable  impact  on  the  ability  to  estimate  the  elevation. 
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Appendix  A.  Nomenclature 


A.l  Symbols 

£  Complex 

*  Complex  conjugate 

=  Defined  as 

G  Element  of 

(•)  Ensemble  average  or  expectation  value 
f  Hermitian  adjoint  (complex  conjugate  transpose) 

In  Identity  matrix,  n  x  n 
Q  Imaginary 

~  On  the  order  of 

5ft  Real 

T  Transpose 

A.2  Acronyms 


AOA 

Angle  of  arrival 

AOB 

Angle  of  bearing 

CRLB 

Cramer-Rao  lower  bound 

FI 

Fisher  information 

LOS 

Line  of  sight 

MCF 

Mutual  coherence  function 

MLE 

Maximum  likelihood  estimator 

MSE 

Mean-square  error 

PDF 

Probability  density  function 

PW 

Plane  wave 

SAA 

Small  angle  approximation 

SNR 

Signal-to-noise  ratio 
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SW  Spherical  wave 

2D  Two  dimensions  or  two-dimensional 

3D  Three  dimensions  or  three-dimensional 
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